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*3.  ABSTRACT 


A  simple  algorithm  this  describes  here  for  removing  nonphysical  self  forces  from 
two  popular  electromagnetic  plasma  simulation  models.  This  algorithm  also  has  two 
additional  features;  it  considerably  reduces  short-wavelength  noise  and  unwanted  numer¬ 
ical  fluctuations,  and  permits  faster  integration  of  the  particle  orbit  equations  by  roughly 
a  factor  of  two. 
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This  note  describes  a  simple  algorithm  for  removing  nonphysical 
self  forces  from  two  popular  electromagnetic  plasma  simulation  models. 
This  algorithm  also  has  two  additional  features ;  it  considerably  reduces 
short-wavelength  noise  and  unwanted  numerical  fluctuations,  and  permits 
faster  integration  of  the  particle  orbit  equations  by  roughly  a  factor 

of  two.  It, is  currently  being  included  in  the  CLYRAB  code  [l ,2], 

V  . 

There  are  three  major  numerical  models  in  which  electromagnetic 
radiation  fields  are  self -consistently  coupled  to  the  Lorentz-liewton 
equations  for  the  charged-particle  motion.  The  first  of  these,  proposed 
by  Buneman  [3],  includes  an  explicit  leapfrog  algorithm  for  solving  the 
Maxwell  Equations  and  a  special  charge-current  algorithm,  based  on  the 
NGP  (Nearest  Grid  Point)  particle  interpolation,  which  ensures  that  the 
continuity  equation  for  charge  and  current  is  automatically  satisfied 
at  each  timestep  in  finite-difference  form.  The  second  of  these  three 
algorithms  [l]  uses  essentially  the  same  field  treatment  but  has  a  more 
flexible  treatment  of  current  accum.lation  in  which,  however,  a  Poisson 
Equation  must  be  solved.  The  third  algorithm  [l]  differs  from  the 
first  two  in  the  treatment  of  the  electromagnetic  fields.  Rather  than 
a  finite  difference  approximation  to  the  Maxwell  Equations ,  the  algorithm 
solves  the  Fourier  transform  of  the  equations  in  k-space.  This  third 
algorithm  has  nut  been  implemented  in  multi dimens ions  but  a  fore¬ 
runner  method  in  one  dimension  was  devised  by  Langdon  and  Dawson  [t j. 

Many  variations  of  the  first  two  algorithms  are  possible  and  several 
have  been  coded  and  used  [5 ,6].  Two  of  these  are  the  !  .-se-Nielson 
algorithms  A  and  B.  Their  algorithm  A  is  an  important  generalization 
of  Buneman' s  algorithm  to  the  PIC  (particle-in-cell)  linear 
interpolation  [8,9].  This  algorithm  also  satisfies  Poisson's 


Equation  automatically  at  every  timestep.  The  Morse-Nielson  algorithm 
B,  as  pointed  out  by  Langdon  [7],  is  essentially  equivalent  to 
the  Boris  algorithm.  Thus  our  comments  apply  to  the  Boris  algorithm 
used  in  CYLRAD,  to  both  Morse-Nielson  algorithms  and  to  the  Sinz  [5] 
algorithm. 

These  electromagnetic  algorithms  treat  particles  moving  on  staggered, 
interlaced  grids  of  variables  as  shown  in  Fig.  1  for  two  dimensions.  Each 
field  variable  is  so  positioned  that  the  time-dependent  Maxwell  Equations 
reduce  to  an  extremely  simple,  finite-difference  form  for  advancing 
£  and  £  which  is  reversible  and  second-order  accurate.  The  generalization 
of  this  2D  field-variable  arrangement  to  3®  is  straight  forward  and  the 
specialization  to  ID  is  trivial.  As  can  be  seen,  the  finite-difference 
form  of  the  divergence  equation  which  should  be  satisfied  at  each  cycle 
is 


=  Unp(i,j).  (1) 

The  electrostatic  fields  which  result  from  placing  a  particle  at  rest 
>>n  the  (i,j)  grid  point  are  different  from  those  fcuna  in  an  electro¬ 
static  code  because  the  x  and  y  electric  field  grids  are  displaced,  as 
seen  in  Fig.  1,  by  Ox/2  and  6y/2  respectively  from  the  positions  they 
would  occupy  in  a  standard  electrostatic  code  [8,9].  To  simplify  the 
analysis,  the  electrostatic  and  the  electromagnetic  grids  are  compared 
in  Fig.  2  for  a  ID  case  where  a  test  particle  of  unit  charge  is  at- 


2 


position  a  off  a  grid  point  (0  <  a  <  6x/2  =  l/2). 

Consider  first  the  electromagnetic  grid  as  used  in  current  simulation 
codes  (Fig.  2).  Hie  linearly  inter j>olated  charge  density  at  grid  points 
0  and  1  are  p(0)  =  1-a  and  p(l)  =  a.  If  +  E^y^*  —  E3/2  are  c*e*>i-ne<*  to 
be  the  electric  field  at  x  *  +  l/2  and  x  *  3/2  respectively  due  to  a 

unit  charge  at  x  —  0,  then  summing  the  contributions  from  the  two  grid 
points  one  gets 

\<-l/2)  *  h/2  -  “  Sl/2> 

Ex(l/2)=  (l-a)  El/2  -  a  ki/2,  (2) 

EX(3/2)  *  <1-a)  h/2  *  °  *1/2- 

Since  the  particle  lies  between  E.^l/2)  and  £^(-1/2),  the  linearly  inter¬ 
polated  x  electric  field,  as  would  be  found  by  an  electromag  etic  simulation 
code,  is 


E  (u)  *  (1/2+  a)  E  (l/2 )  +  (1/2  -  u)  E  (l/2 ) 


*  **  [El/2-  ^  b/2 


t  0. 


b) 


Equation  (5)  shows  that  the  self  electrostatic  field  of  a  single  simulation 
parti ;le  is  non-zer Thus  all  sorts  of  spurious  effects  can  result. 

Figure  3  shows  a  ID  plot  of  the  equivalent  potential  a  particle  woxild  see 
due  to  its  self-force.  We  have  carried  out  tests  on  an  electromagnetic 
code  and  the  oscillations  of  a  particle  in  this  self  field  have  been 
observed.  The  preceeding  analysis  has  been  generalized  to  electrostatic 
and  magnetostatic  self-forces  in  two  and  three  dimensions.  In  every 
case  the  results  are  the  same;  spurious  electrostatic  and  magnetostat'C 


forces  are  found  when  the  charges  do  not  exactly  lie  on  grid  lines.  We 
know  that  there  are  no  self  forces  in  the  usual  electrostatic  codes  [6,9] 
and  the  electrostatic-code  field  in  Fig.  2  are  clearly  closely  related 
to  the  electromagnetic-code  fields.  Therefore,  it  should  be  a  simple 
matter  to  eliminate  the  electrostatic  self-forces  from  the  electromagnet!' 
code.  From  Fig.  2  we  have 


E^l/2)  =  <t(l)  -  $(C) 
6x 


(M 


and 


E^O)  =  tt(l)  -  0(-l) 

5x 


b) 


the  fields  are  therefore  related  by 


£*(0)  a  1/2  [£^(1/2)  +  E^-l/2)] 


(6) 


If  we  linearly  interpolate  to  the  particle  position  using  the  averaged 
Helds  of  Eq.  (6)  and  Eq.  (2),  we  get 


E^a)  =  (1-a)  Ex(0)  +  a  Ey(l) 


(7; 


=  0. 


This  is  the  desired  result  of  zero  self -force.  When  generalized  to  two 
and  three  dimensions,  the  electrostatic  and  magnetostatic  self  forces 
are  found  to  be  zero.  Furthermore*,  since  the  argument  is  basically  ~ne 
of  symmetry  rather  than  being  based  on  any  particular  force  law,  the 
determination  of  $>(i)  from  p(i)  admits  all  finite-sized  particle  algorithms 


...  . 


as  well  as  the  3 ,5,  and  7  point  Poisson  operators  found  in  one,  trio, 
and  three  dimensions. 

The  averaging  alvoriti  _  developed  here  for  fully  electromagnetic 
simulations  in  two  ~  i/n  ..nsions  is  an  obvious  extension  of  Eq.(6)to  the 
field-veriabi,  e  layout  shown  in  Fig.  1.  Die  three  components  of*  particle 
current  density,  J^(i,j),  J^(i, j )  and  J^(i, j )  are  all  found  by  linear 
interpolation  onto  the( i, j )grid,  exactly  as  p(i,j ) .  The  current  densities 
for  use  or  the  interlace^  field  grids  are  then  computed  as  follows: 

J  t  l-l  -  J)  5  1L  [J?(i,j)  + 

JyVi,j+l/2)  =  1/2  [j£(i,.i)  +  j£(i,o+l)3>  (o) 

The  six  field  components,  F  ,  E  ,  E  ,  B  ,  3  ,  and  B  are  then  all  in- 

x  y  z  x  y  z 

tegrated  exactly  as  prescribed  by  the  staggered-leapfrug  algorithm. 

These  field  components  are.  however,  averaged  bacx  to  the  particle  grid 
before  being  used  to  advance  the  particle  equations  of  motion.  Thus 

liJCij)  =  1/2  rEx(i+l/2i,.i)  +  Ex(i-l/2,jj3 

=  1/2  [Ey(i,j+l/2)  +  Ey(i,j-l/2)] 

E^(i,j)  E  Ez(i,j), 

)  =  1/2  [Bx(i,o+l/2)  +  Bx(i,j-l/2)],  {' ! 

Bj(i,j)  =  1/2  [Bz(i+l/2,j)  +  Bz(i-l/2,j)3, 

B2(i,j)  =  lA  [Bz(i+l/2,itl/2)  +  Bz(i+l/2,j-l/2) 

+  Bz(i-l/2  ,j+l/2)  +  Bz(i-l/2„j-l/2)3. 
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The  original  fields  JS  and  B  are  retained  unchanged,  however,  for  use  in 
the  next  cycle  to  advance  Maxwell's  Equations.  They  are  not  to  be  computed 
as  averages  of  the  particle  fields  jgPand  Bp. 

The  immediate  consequence  of  adding  the  averaging  stages  given  by 
Eqs.  (8)  and  (9)  to  standard  fully  electromagnetic  codes  is  the  removal 
of  nonphysical  electrostatic  and  magnetostatic  self -forces.  There  are 
two  other  favorable  and  important  consequences.  First,  the  averages 
required  are  smoothing  operation?;  therefore,  spurious  numerical  Cherenkov 
radiation  and  brenmstrahlung,  arising  mostly  at  short  wavelengths  should 
be  strongly  suppressed.  Second,  a  sizeable  simplification  results  since 
all  particle  quantities  are  now  defined  on  a  single  grid,  only  one  set 
of  bilinear  weight  coefficients  need  be  found,  rather  than  four,  and  thus 
it  is  expected  that  optimized  particle  integration  can  be  speeded  up  by 
at  least  a  factor  of  two. 
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